/* *****************************************************************************
*
* Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
* All Rights Reserved.
*
* This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
* the contents of this file may not be disclosed to third parties, copied or
* duplicated in any form, in whole or in part, without the prior written
* permission of Silicon Graphics, Inc.
*
* RESTRICTED RIGHTS LEGEND:
* Use, duplication or disclosure by the Government is subject to restrictions
* as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
* and Computer Software clause at DFARS 252.227-7013, and/or in similar or
* successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
* rights reserved under the Copyright Laws of the United States.
*
***************************************************************************** */
      subroutine add_circle( nlin, ncol, mat, lda)
      real mat(lda,ncol)

	ior = nlin / 2
	jor = ncol / 2

	r = MIN(ior,jor)/9.
	r2 = r * r

c$doacross local(i,j,dx,dy,d2,z)
	do j = 1, ncol
	    do i = 1, nlin
		dy = j - jor
		dx = i - ior
		d2 = dx*dx + dy*dy
		z = (.01 + d2) / r2
		mat(i,j) = Sin(z) / z
c		mat(i,j) = Sin(z) * Cos(dx/r) / z
	    end do
	enddo

      return
      end


      subroutine inifil(a,lda,nlin,ncol)
      real a(lda,1)
c************************************************
c						*
c   create a wave propagation filter		*
c						*
c   The Fourier Transform of the 2D wave	*
c   equation link the "time" frequency W, with	*
c   the "space" frequencies KX and KY		*
c						*
c	W^2 == c^2 * (KX^2 + KY^2), c velocity	*
c						*
c   The wave propagation between time 0 and T 	*
c   simulated is given by the phase shift ope-	*
c   rator   Exp( - i * W * T)			*
c   which translate in :			*
c	Exp( - i * c * Sqrt(KX^2+KY^2) * T	*
c   if we consider the forward propagation,	*
c   and its conjugate for the backward		*
c						*
c************************************************
c						*
c   Here we choose :				*
c	c = 1					*
c	T = 1/ ncol				*
c						*
c************************************************
	i_lin_even= mod(nlin+1,2)
	i_col_even= mod(ncol+1,2)
	T = 1 / Float(MIN(ncol,nlin))
c
c   first line  ( KX == 0 )
c
	KX = 0.
	a(1,1) = 1.0
	do KY = 1, (ncol-1)/2
	    y = Float(KY)
	    W = Sqrt( y*y )
	    a(1,2*KY) = cos(W*T)
	    a(1,2*KY+1) = - sin(W*T)
	end do
	if( i_col_even) then
	    y = Float(ncol/2)
	    W = Sqrt( y*y )
	    a(1,ncol) = cos(W*T)
	endif
c
c   General Values
c
	do KY = 0, (ncol-1)
	    y = Float(KY)
	    do KX = 1, (nlin-1)/2
		x = Float(KX)
		W = Sqrt(x*x + y*y)
		a(2*KX,KY+1) = cos(W*T)
		a(2*KX+1,KY+1) = - sin(W*T)
	    end do
	end do
c
c   first line  ( KX == nlin/2 )
c
      if( i_lin_even ) then
	KX = (nlin/2)
	x = Float(KX)
	a(ncol,1) = 1.0
	do KY = 1, (ncol-1)/2
	    y = Float(KY)
	    W = Sqrt( y*y + x*x)
	    a(ncol,2*KY) = cos(W*T)
	    a(ncol,2*KY+1) = - sin(W*T)
	end do
	if( i_col_even) then
	    y = Float(ncol/2)
	    W = Sqrt( y*y + x*x)
	    a(ncol,ncol) = cos(W*T)
	endif
      endif

      return
      end
      
      subroutine reverse(a,lda,nlin,ncol)
      real a(lda,1)
	i_lin_even= mod(nlin+1,2)
c
c   first line  
c
	do j = 3, (ncol-1), 2
	    a(1,j) = - a(1,j)
	end do
c
c   General Values
c
	do j = 1, ncol
	    do i = 3, (nlin -1), 2
		a(i,j) = - a(i,j)
	    end do
	enddo
c
c   last line  
c
	if( i_lin_even) then
	    do j = 3, (ncol-1), 2
		a(nlin,j) = - a(nlin,j)
	    end do
	end if
      return
      end
      
      subroutine matnoise(a,lda,nlin,ncol)
      real a(lda,1)
c
      init = 1325

      do 30 j = 1,ncol
         do 20 i = 1,nlin
            init = mod(3125*init,65536)
            a(i,j) = a(i,j) + (init - 32768.0)/16384.0
   20    continue
   30    continue
      return
      end
      
      subroutine matpower( nlin,ncol,mat,power,lda)
      real mat(lda,*), power(lda,*)
	real re, im

	    il = 1
	    power(il,1) = mat(il,1) * mat(il,1)
	    do j = 2, (ncol-1)/2
		re = mat(il,2*j-1)
		im = mat(il,2*j)
		power(il,j) = re * re + im * im
		power(il,ncol-j+1) = power(il,j)
	    end do
	    if( mod(ncol,2) .eq. 2) 
     $		power(il,ncol) = mat(il,ncol) * mat(il,ncol)
c$doacross local(i,j,re,im)
	do j = 1, ncol
	    do i = 2, (nlin-1)/2
		re = mat(2*i-2,j)
		im = mat(2*i-1,j)
		power(i,j) = re * re + im * im
		power(nlin-i+1,j) = power(i,j)
	    end do
	end do
	    il = nlin
	    power(il,1) = mat(il,1) * mat(il,1)
	    do j = 2, (ncol-1)/2
		re = mat(il,2*j-1)
		im = mat(il,2*j)
		power(il,j) = re * re + im * im
		power(il,ncol-j+1) = power(il,j)
	    end do
	    if( mod(ncol,2) .eq. 2) 
     $		power(il,ncol) = mat(il,ncol) * mat(il,ncol)
	if ( mod(nlin ,2) .eq. 0) then
	end if
      return
      end

      subroutine matlog( nlin, ncol, power, lda)
      real power(lda,*)
	real t
c$doacross local(i,j,t)
	do j = 1, ncol
	    do i = 1, nlin
		t = power(i,j) + 1
		power(i,j) = LOG(t)
	    end do
	end do
      return
      end
	    
      subroutine matcolor( nlin,ncol,mat,color,lda)
      real mat(lda,*)
      integer color(nlin,*)
	real t, a, coeff,xx,yy
	integer red, green, blue

	xx = 0.0
	yy = 0.0
	do j = 1, ncol
	    do i = 1, nlin
		a = mat(i,j)
		if ( a .lt. xx) then
		    xx = a
		else if(  a .gt. yy) then
		    yy   = a
		end if
	    end do
	end do
	coeff = MAX(ABS(xx),ABS(yy))

c	print *, ' MIN =', xx, '   MAX =', yy

	if ( coeff .gt. 0.) coeff = .99 * 255. / coeff

c$doacross local(t, a, red,blue,green,i,j)
	do j = 1, ncol
	    do i = 1, nlin
		t= mat(i,j)
		a = ABS(t)
		red = coeff* (a - t)
		if( red .gt. 255) red = 255
		green = coeff* a
		blue = coeff* (a + t)
		if( blue .gt. 255) blue = 255
		color(i,j) = (red) + (256*green) + (256*256*blue)
	    end do
	end do
      return
      end
	    
